#include "libmesh/quadrature.h"
#include "Assembly.h"
#include "DensityJumpIndicator.h"
#include "CFDProblem.h"

template<>
InputParameters validParams<DensityJumpIndicator>()
{
  InputParameters params = validParams<FluxJumpIndicator>();
  return params;
}

DensityJumpIndicator::DensityJumpIndicator(const InputParameters &parameters) :
	FluxJumpIndicator(parameters)
{
}

Real DensityJumpIndicator::computeQpIntegral()
{
	for (size_t i = 0; i < _uh.size(); ++i)
	{
		_cfd_data.uh[i] = (*_uh[i])[_qp];
		_cfd_data_neighbor.uh[i] = (*_uh_neighbor[i])[_qp];
	}
	_cfd_data.reinitInviscous();
	_cfd_data_neighbor.reinitInviscous();
	Real vn = (_cfd_data.vel+_cfd_data_neighbor.vel)*_normals[_qp] + _cfd_data.c+_cfd_data_neighbor.c;

	return fabs((_cfd_data.r - _cfd_data_neighbor.r)/(_cfd_data.r + _cfd_data_neighbor.r)*vn);
}

